#include "mex.h"
#include <math.h>
#include <assert.h>

void Win2Ws(double *Ws, double *Win, double *P, double *T);

void mexFunction(
	int nlhs, mxArray *plhs[],
	int nrhs, const mxArray *prhs[])
// Matlab call: Ws = Win2Ws(Win,P,T)
{
	double  *Ws;
	double	*Win;
	double	*P;
	double	*T;
	
	if(nrhs!=3)
		mexErrMsgTxt("Three inputs required");
	if(!mxIsDouble(prhs[0]) || mxIsComplex(prhs[0]) ||
	   !mxIsDouble(prhs[1]) || mxIsComplex(prhs[1]) ||
	   !mxIsDouble(prhs[2]) || mxIsComplex(prhs[2]))
		mexErrMsgTxt("Inputs must be a scalar");
	if(nlhs>1)
		mexErrMsgTxt("Only one output is allowed");
		
	Win = mxGetPr(prhs[0]);
	P   = mxGetPr(prhs[1]);
	T   = mxGetPr(prhs[2]);
	
	plhs[0] = mxCreateDoubleMatrix(1,1,mxREAL);
	Ws  = mxGetPr(plhs[0]);
	
	Win2Ws(Ws,Win,P,T);
}

void Win2Ws(double *Ws,double *Win, double *P, double *T)
{
	double Xin;
	double Xs;
	double Win_temp;
	double Win_max;
	
	Xin = 0.001;
	Win_temp = 5.7423e4*(1.6e4*exp(-18500/(*T))*((*P)*Xin))/(1+0.19*sqrt((1-Xin)*(*P))+0.70*(*P)*Xin);
	Win_max = 5.7423e4*(1.6e4*exp(-18500/(*T))*((*P)*1.0))/(1+0.19*sqrt((1.0-1.0)*(*P))+0.70*(*P)*1.0);
	assert((*Win)<=Win_max);
	while(Win_temp < (*Win))
	{
		Xin += 0.001;
		Win_temp = 5.7423e4*(1.6e4*exp(-18500/(*T))*((*P)*Xin))/(1+0.19*sqrt((1-Xin)*(*P))+0.70*(*P)*Xin);
	}
	if(Xin > 0.4377)
	{
		Xin = 0.4377;
	}
	Xs = Xin/0.4377;
	(*Ws) = 5.7423e4*(1.6e4*exp(-18500/(*T))*((*P)*Xs))/(1+0.19*sqrt((1-Xs)*(*P))+0.70*(*P)*Xs);
}
